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ABSTRACT 

We combine the the stellar spectral synthesis code Starburst99, the nebular modelling code MAP- 
PINGS III and a 1-D dynamical evolution model of Hii regions around massive clusters of young 
stars to generate improved models of the spectral energy distribution (SED) of starburst galaxies. 
We introduce a compactness parameter, C, which characterizes the specific intensity of the radiation 
field at ionization fronts in Hii regions, and which controls the shape of the far-IR dust re-emission, 
often referred to loosely as the dust "temperature" . We also investigate the effect of metallicity on 
the overall SED and in particular, on the strength of the PAH features. We provide templates for 
the mean emission produced by the young compact Hii regions, the older (10 — 100 Myr) stars and 
for the wavelength-dependent attenuation produced by a foreground screen of the dust used in our 
model. We demonstrate that these components may be combined to produce a excellent fit to the 
observed SEDs of star formation dominated galaxies which are often used as templates (Arp 220 and 
NGC 6240). This fit extends from the Lyman Limit to wavelengths of about one mm. The methods 
presented in both this paper and in the previous papers of this series allow the extraction of the 
physical parameters of the starburst region (star formation rates, star formation rate history, mean 
cluster mass, metallicity, dust attenuation and pressure) from the analysis of the pan-spectral SED. 
Subject headings: ISM: dust — extinction, H ii — galaxies:general,star formation rates, starburst- 
infrared : galaxies — ultraviolet : galaxies 
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1. INTRODUCTION 

By definition, the bolometric luminosity of a starburst 
galaxy is dominated by the young stars it contains. Thus 
regardless of how much or how little of this luminosity 
is reprocessed through the dusty interstellar medium ei- 
ther through thermal emission in the IR of dust grains, 
through fluorescent processes, or through heating and 
re-emission in an ionized medium, the pan-spectral SED 
encodes information about what the star formation rate 
currently is, and what it has been in the recent past. The 
first objective of pan-spectral SED modelling is therefore 
to be able to reliably infer star formation rates in galax- 
ies and to provide likely error estimates using observa- 
tional data sets which may in practice be restricted to 
only certain emission lines or spectral features. In prin- 
ciple, almost any part of the SED of a starburst can be 
used as a star formation indicator, provided that the ap- 
propriate bolometric correction to the absolute luminos- 
ity can be made, and observational issues are accounted 
for. In practice, each wavelength regime has a different 
level of sensitivity to the ongoing star formation which 
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is dependent upon these bolometric corrections, with hy- 
drogen emission lines and the IR part of the SED being 
the most robust indicators of the current star formation 
rate (SFR). These bolometric corrections critically de- 
pend on the foreground dust absorption (more properly 
called dust attenuation), and the geometry of the em- 
bedded dusty molecular clouds, with respect to both the 
ionizing stars and the older stellar population. The ac- 
curate determination of such bolometric corrections is a 
major motivation of our theoretical work of pan-spectral 
SED modeUing. 

In order to correctly model the SEDs of starburst 
galaxies, we first need to understand how the form of the 
SED is controlled by the interstellar physics and the ge- 
ometry of the stars with respect to the gas. Once these 
are understood, we can then use our theoretical mod- 
els to attain the objectives of our second motivation for 
such SED modelling; to gain insights into the physical 
parameters of starburst galaxies. In particular, we can 
hope to quantify the stellar populations, the atomic and 
molecular gas content, the star and gas-phase metallici- 
ties, physical parameters of their interstellar media such 
as the pressure or mean density, and the nature of the 
interstellar dust, both its composition and spatial distri- 
bution. The physical parameters so derived on homo- 
geneous samples of objects can then help develop our 
insight into the physical processes which control them. 

The dust grain temperature distribution, and therefore 
the shape and peak of the far-IR feature, depends crit- 
ically on the geometrical relationship between the dust 
grains and the stellar heating sources assumed in the 
model. Models with warmer far-IR colors will have a 
more compact disposition of gas with respect to the stars. 
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The difficulty here is that, in any simple starburst model, 
these geometrical relationships are not determined a pri- 
ori. 

In the semi-empirical modelling of Dale and his collab- 
orators f Pale et al.l[20?)ll : iDale & Helod[2?)M ). the SEDs 
of both disk and starburst galaxies were suggested to 
form a one-parameter family in terms of dust tempera- 
ture. This suggested th at starburst galaxies ha v e hot - 
ter dust temperatures. iLagache. Dole fc Pugetl ()2003D 
(again empirically) have suggested that the absolute lu- 
minosity controls the form of the SED. Both of these as- 
sertions may be true to some extent, since IR luminous 
galaxies generally have greater rates of star formation 
than normal galaxies.^ 

The French group (!Galli ano et al.|[2003f ) take the sim- 
plest approach of approximating the starburst by a 
spherical H II region and clumpy dust shell around 
the central star for ming region. A more ad v anced ap- 
proach is used by ISiebenmorgen fc Kriigel ()2007l ) to 
model starbursts and ULIRGs. While assuming a spher- 
ical geometry for the radiative transfer, they also in- 
clude the effect of the hot dust around young OB-stars 
as well as the diffuse ISM dust surrounding the star- 
burst and older stellar population. The hot dust compo- 
nent is important as it can dominate the mid-IR emission 
pCriig el fc Siebcnmorgen 1993). 

Associated with the latter models is the approach by 
lEfstathiou et al.l (|200Cl( ). Likewise concentrating on star- 
burst galaxies, they modelled the starburst as a group of 
stars surrounded by thick molecular clouds in a manner 
similar to our approach, and in addition, also included 
a similar, simple description for the evolution of the dis- 
tance of the molecular clouds to the illuminating stars. 
With these models they could explain the observed IRAS 
distributions and reproduce several ISO observations. 

One of the more sophisticated approaches is taken 
in the GRASIL code by the Pado va-Trieste group 
(jSilva et al.lll998t iGranato et al.ll2000D . Their starburst 
model uses a spherical geometry with King profiles, 
and they allow for the formation of clusters of stars in 
molecular complexes, and their subsequent escape from 
these regions. This group has sin ce incorporated gas 
physics by use of the CLOUDY code (jFerland et al.lll998l : 
I Abel et ani2005f ) to provide emissi on line diagnostics a s 
well as dust continuum diagnostics (jPanuzzo et al.|[2003l ). 

A similar ly adva nced approach was used by 
iPiovan et all ()2006al lbl) . who also assume a spheri- 
cal geometry with King profiles and young stars in 
molecular complexes. 

In the conceptually sophisticated models of 
iTagaki et all C2003a,b), ' a mass-radius relationship for 
the star formation region of r^/kpc = 0(M/1O^M0)-'^/^ 
is adopted along with a stellar density distribution 
given by a generalized King profile. The parameter 8 
is a compactness parameter which expresses the degree 
of matter concentration, and is related to the optical 
depth of the dust through which the starburst region is 
seen. For a sample of ultra-luminous starbursts, they 
find that, while most conform to a constant surface 
brightness of order lO^^LQkpc"^, there are a few objects 
with surface brightnesses roughly ten times larger than 
this, which they ascribe to post-merger systems. In this 
paper, we will adopt a derivative version of this concept 



of a compactness parameter as the factor which provides 
the main control on the shape of the far-IR bump. 

All the fully theoretical (as opposed to semi-empirical) 
methods used by other groups involve the calculation of 
essentially a single spherical radiation transfer problem. 
Unfortunately, real starburst galaxies have many sepa- 
rate clusters of many different ages distributed in a com- 
plex spatial distribution. However, gas column densities 
and pressures can be extremely high, and as a conse- 
quence the size scale of individual H II region complexes 
can be extremely small in comparison to the overall scale 
of the starburst. It is only when the star forming com- 
plexes join up to produce large-scale collective phenom- 
ena as in, for example, the outflow in M82, that we need 
to go to a full 3-D radiative transfer model covering the 
whole galaxy. 

In the earlier papers in this series, we take advantage 
of the localized radiative transfer approximation to con- 
struct our pan-spectral SEDs. Instead of treating the 
starburst as a single H II region complex covering the 
whole starburst region, we split the starburst up into 
many individual H II regions, each ionized by the UV 
photons of the clusters within them, and each evolving 
in radius and internal pressure according to the mechani- 
cal energy input of the exciting stars through their stellar 
winds and supernova explosions. The global SED is then 
the sum of the SEDs produced by each of these H II 
regions and their surrounding photo-dissociation regions 
(PDRs) integrated over all cluster ages. Since the radia- 
tive transfer problem in each H II region is fully treated, 
in this approach we stand a better chance of capturing 
the full range of physical conditions encountered in a 
starburst region. This collective approach is what can be 
found in many of the sophistic ated modelhng codes such 
as that of lPiovan etHI (|2006af) and GRASIL (jSilva et all 
1998). In GRASIL for example, multiple core molecular 
cloud systems can be included with differing parameters 
for each, such as mass and optical depth. 

In paper I (|Dopita et al.ll2005al) . we investigated the 
role that pressure alone plays in changing the compact- 
ness of the H II regions within the starburst and hence 
in the controlling the shape of the far-IR dust emission 
bump. However, a defect in these models is that they 
were only run with a single value of mean cluster mass 
(Mci). A change in cluster mass directly affects the spe- 
cific intensity of the radiation field in the H II region, and 
this will in turn chang e the shape of the fa r-IR bump. 

In papers H & HI (jPoDita et all I20n6al lbl). we intro- 
duced the TZ = (Mci) /Pq parameter which controls the 
absolute value of the ionization parameter in the H II re- 
gion and its time evolution. The ionization parameter is 
defined as the ratio of the ionizing photon density to the 
particle density in the H II region; lA = L\jy /AnR^^^nc, 
where Lyv is the flux of ionizing UV photons produced 
by the central cluster, i?Hii is the mean radius of the 
H II region with particle density n, and c is the speed of 
light. All models having a given value of TZ and metal- 
licity Z will show the same run of ionization parameter 
as a function of cluster age and will therefore produce 
identical line ratios at any given age. 

The shape of the dust feature or "bump" in the IR 
in starbursts is controlled by the distribution of dust 
temperatures in the starburst galaxy. Within a given 
H II region, this distribution of temperatures is con- 
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trolled by the specific photon density, meaning that 
the mean dust temperature of any individual grain is 
(Tgi) = /(Luv/^Hii)- Thus denser and more compact 
H II regions will produce hotter grain temperature dis- 
tributions. In this paper, we introduce, by analogy with 
the TZ parameter, a "Compactness Parameter", C. All 
models having a given value of C and metallicity Z will 
show the same run of grain temperature distribution as 
a function of cluster age and will therefore produce iden- 
tical far-IR dust re-emission bumps at any given age. 

In the following sections of this paper, we will discuss 
the details of our modeling procedure, insofar as this is 
different from that used in earlier papers in this series, in- 
troduce the Compactness Parameter and show that this 
does indeed serve to characterize the far-IR bump for H II 
regions in the age range 0.5 — 10 Myr. We also investigate 
the effects of varying metallicity on the form of the far- 
IR bump, provide templates for the mean SED of com- 
pact H II regions, and of the older (10 — 100 Myr) stars, 
and for the attenuation produced by a dusty fractal fore- 
ground screen. Finally, we show how these components 
can be combined to produce excellent fits to frequently- 
used star burst templates such as Arp 220 or NGC 6240. 

2. MODELS 

The Starburst models calculated here follow the gen- 
eral form described in the p revious papers of this series 
(|DoDita et al.ll2005al [2006al in): hereinafter SEDl, SED2 
and SED3, respectively. However, apart from being up- 
dated to use the latest versions of the modelling codes 
Starburs t99 (Lcithcrcr et al. 1999) and MAPPINGS ill 
([Groves! 12004 ) . the models incorporate a number of 
changes or improvements that we here describe in greater 
detail. To make this paper self-contained, we briefly re- 
capitulate on the techniques used in the earlier papers of 
this series. 

2.1. Stars and Stellar Clusters 

We have used the latest version (2006) of Starburst99^ 
to compute the spectral energy distribution of clusters 
of stars of any given age. A detailed description of 
latest stellar atmospheres and s tellar evolut i on ph ysics 
used within the code are given in lSmith et al.l ()2002l ) and 
IVazauez fc Leithereij ([2005). 

In our Starburst99 models, we t ake an in stantaneous 
burst of Mci = IO^Mq, having a iKroupal f2002 ) bro- 
ken power-law IMF between 0.1 and 120 Mq. Within 
the code we use the standard combination of the 
Geneva and Padova track s for the stellar evolution 
(jVazquez fc Le ithcrcr 200^ , and these determine the to- 
tal mechanical luminosity, Lmcch- We use the theoretical 
"high-mass-loss" tracks for the treatment of the stellar 
wind. The supernova cutoff mass is 8 Mq. However, 
since the H II region evolution is run up to only 10 Myr, 
the exact choice of this cutoff mass is unimportant to the 
modelling. 

We output the stellar spectra at 0.01 Myr, 0.5 Myrs 
and then every 0.5 Myrs after that up to an age of 10 
Myrs, by which effectively all of the ionizing photons have 
been emitted (see SED2). Note that the resolution of the 
starburst model is actually higher at 0.1 Myrs, which 
provides the fine-gridding needed to accurately track the 



mechanical energy input used in the computation of the 
evolution of the H II regions. For completeness, this 
model was also extended to 1 Gyr, with spectra com- 
puted at longer intervals for older stellar templates. 

To determine the parameters for clusters of any given 
mass, we assume a simple scaling with cluster mass. Such 
a scaling should hold in starbursting regions where many 
clusters are forming, as here stochastic effects are gen- 
erally small, and on average the IMF is well sampled 
throughout the mass range. However, as a number of au- 
thors have pointed out, sever al of whom are refere nced in 
SED2, and most recently bv lWeidner fc Kroupal (|2006), 
such an assumption will not hold for low cluster masses. 

2.2. Hu Region Evolution 

The H II regions are treated as 1-D mass-loss bubbles 
driven by th e mechanical energ y input of their stars and 
supcrnovae (jCastor et al.|[l975r ). Their equation of mo- 
tion is given by (SEDl): 
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where the time dependent mechanical luminosity, 
Lmcchit), is determined from the Starburst99 output. 

The pressure in the H ii region with radius i?, expand- 
ing in an ISM with density po is then determined as, 
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where uhii is the density in the ionized gas, with elec- 
tron temperature Tg . This equation is deri ved from the 
Oev fc Clark"a[l997l [l998h version of the ICastor et al.l 



mass-loss bubble formulae with the assumption 



that the H ii region has the same pressure as the shocked 
stellar wind, and is confined to a thin shell around the 
periphery of the wind-blown bubble. The ionizing flux at 
the inner boundary of the Hii region is then Ljjv /4:TtR'^ 
and the density in the ionized region uhii is given by 
equation 12. 2( from which the ionization parameter of the 
H II region, U, can be derived. 

2.3. Nebular Abundances and Depletion Factors 

The abundance set and depletion factors used in these 
models are unchanged from those presented in SED3 and 
are given in table [Tl Th e nebular abundance set fol- 
lows |As2lunTietIal] ()2005l ). As noted previously, the gas 
phase "Solar" abundance in the models is somewhat off- 
set from the "Solar" abundance set used in Starburst99. 
While this has been shown to have no significant effect 
on the models, it does result in a small inconsistency in 
the models. 

The exploration of metallicity effects within the mod- 
els presented here is limited to those metallicities com- 
puted in Starburst99; 0.05^0, 0.2Zq, 0.4^0, I.OZq, and 
2.0Zq. As in SED3, we simply scale the abundances with 
metallicity, with the following exceptions. For helium we 
use the empirical relationship to include the primordial 
component as well as that from nucleosynthesis: 



http:/ /www. stsci.edu/sciencc/starburst99/ 



He/H ^ 0.0737 + 0.024(Z/Zc 



(3) 
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The elements carbon and nitrogen are both observed to 
have a primary nucleosynthetic component and a sec- 
ondary nucleosynthetic component at higher metalhci- 
ties. Thus for these elements we adopt the following em- 
pirical relationships; 

C/H = 6.0 X 1O-^(Z/Z0) + 2.0 X X^^-'^^ZIZ^f (4) 
N/H=l.l X lO"^(Z/Z0) + 4.9 X 10"'^(Z/Zq)2 (5) 

The depletion factors used in the models are based 
upon the obs erved depletion frac tions in the local inter- 
stellar cloud (' Kimura eraLll200l . These may not repre- 
sent in reality those found in the starburst environments, 
but currently we have no adequate way of estimating 
these. We must also assume that the depletion pattern 
does not vary with metallicity, as there are currently no 
good models or observations for the relationship between 
metallicity and depletion. However, in the starbursting 
environm ents modelled here, such an assumption may be 
adequate (jPraine et al.ll2007l ). With this assumption the 
dust-to-gas ratio is purely a function of metallicity. 

2.4. Photoionization models 

For the component H II regions, we compute both the 
emission and internal absorption of both gas (line plus 
continuum) and dust (continuum including polycyclic 
aromatic hydrocarbons) using the MAPPINGS ill code, 
with the Starburst99 stellar cluster spectra as input. 

Apart from the effects of burst age and metallicity, here 
we explore three other parameters within the models; the 
ISM pressure, Pq, the cloud covering fraction, /pdr, and 
the Compactness Parameter, C, described below, which 
is a function both of Pq and of the mean cluster mass, 
(Mci). 

For the models investigated here we examine five dif- 
ferent thermal gas pressures; Po/k = 10'*, lO'"^, 10^, 10^, 
and 10* K cm~^. These five pressures cover the full range 
expected to be encountered in starbursting galaxies, from 
regions of enhanced star formation in disc galaxies, to the 
high pressure ULIRGs. 

For each parameter set we compute models at 21 star- 
burst ages, covering the timescale 0.01 Myrs to 10 Myrs 
in steps of 0.5 Myrs. 

Finally, for each age we run two models. The first 
model is of the H ii region alone; the region within which 
99% of the hydrogen line emission arises, and within 
which almost all of the ionizing photons are absorbed. 
It is within this region that the hottest dust emission 
arises. The second model corresponds to the Photodisso- 
ciation Region (PDR) surrounding the H II region. This 
region is the transition layer between the H ii region and 
surrounding dense molecular cloud, from which the stel- 
lar cluster is thought to have formed. In the PDR, a 
large fraction of the stellar light is absorbed and most of 
the PAH and dust far-IR emission produced. We fol- 
low the radiative transfer beyond the ionization front 
in the H II region until a total hydrogen column depth 
of A^(HI) = lO'^^ cm^^ is achieved. The column depth 
of 10^^ cm~^ is based upon both observations and upon 
theoretical considerations. The observational data comes 
from measurements of individual molecular cloud s both 
within our galax y (|Larso ^ [l98ll iSolomon et all [19871 : 



iHever et al.l 20011) and in neig hbouring galaxies such as 
M33 ( Rosolowskv et all l2003f) which all give hydrogen 



column densities of ~ lO'^^ cm~^ independent of cloud 
radius. This value is not unexpected, since it indicates 
that all giant molecular clouds are marginally stable 
against gravitational collapse, provided that their virial 
temperatures are a few tens of degrees K. 

These two models, Hii and PDR, are combined 
through our final parameter, /pdr, the starburst cloud 
covering fraction. This parameter is a simplified ver- 
sion of the clearing timescale introd uced in Paper I , 
and discuss ed in other dust models feg'Si lva et al.lll998t 
ICharlot"fc"F all 2000). We introduce this parameter be- 
cause a starbursting system will be a conglomerate of 
bursts of different ages and sizes, unlike a single molec- 
ular cloud around an individual cluster. Thus, while the 
molecular cloud clearing timescale offers a more physi- 
cal picture for an individual cluster, the starburst cloud 
covering fraction better represents what is likely to be 
encountered in a starbursting system. 

The multiple star clusters forming in a starburst will 
be of all possible masses, and, sampled at any instant in 
time, of all possible ages. To account for this we compute 
the luminosity-weighted average of all twenty one ages 
computed between — 10 Myr. In figure [1] we show the 
SEDs of the 21 calculated ages of an individual H ii region 
and a H II region with its PDR, along with the summed 
final average SED for each. These figures clearly show 
the evolution of both the stellar and nebular spectrum 
with age, and reveal the decreasing cluster UV flux and 
cooling dust temperatures as the clusters age and the H li 
bubble expands. 

In all, we have computed a total 300 starburst H II 
region models covering 5 metallicities, 6 values of the C 
parameter (described below), 5 values of the ISM pres- 
sure, and a separate Hii region and H II region plus 
PDR model for each set of parameters (corresponding to 
/pdr ~ and 1, respectively). All models are scaled in 
flux to correspond to a star formation rate of 1 MQyr~^ 
continued over the 10 Myr lifetime of the H II regions. 



2.5. Dust Physics 

The treatment of dust within the photoionization code 
MAPPINGS III was discussed in SEDl. Here, we con- 
centrate only those areas where changes have been made 
to the dust parameters within the code. 

In brief, our dust model consists of 3 components; 
Graphites, Silicates and Polycyclic Aromatic Hydrocar- 
bons. The optical data for each of these comes from 
the series of work by Drainc et al.s (L aor fc Draiii3ll993t 
iLi &: Drain j [2001 ). The IR spectrum arising from dust, 
excluding the effects of PAH emission, is calculated self- 
consistently, including the effects of stochastic heating in 
sn iall grains (al l dust calculations are discussed in detail 
in lGrove^l2004D . The total dust-to-gas ratio within the 
code is set by the fraction of metals depleted from the 
gas onto dust, given in table [TJ 

The heavy elements removed from the gas phase are 
distributed between the two main types of dust, carbona- 
ceous and siliceous, with the carbonaceous dust being 
further divided into graphite and PAHs. The graphite 
and silicate dust is distributed across a grain size distri- 

^ http:/ /www. astro. princeton.edu/~drainc/dust/dust. did. html 
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bution arising from grain destruction processes (SEDl); 

e-(a/a„,„)-^ 

dNia)/da^ka-'-'^-^-^^j^, (6) 

with k defined by the dust-to-gas ratio. 

The minimum grain size of graphite and sihcates are 
cimin =20A and 40A respectively while the maximum 
grain size is the same for both species at Omax = 1600 A. 

2.6. PAHs 

PAHs are treated somewhat separately to the other 
types of dust. They are given a characteristic size and 
an opacity similar to coronene, the best studied catacon- 
densed PAH. We have constructed an empirically-based 
IR emission spectrum described in more detail below. 
The PAH-to-gas ratio is defined through the PAH-to- 
Carbon dust ratio, which is set at 30% in these models. 
While this m ay not be accurate f or all environments and 
metallicities (jPraine et al.l[2007f ). it provides a reason- 
able match to current observations of nearby starform- 
ing galaxies. Differences between the template IR emis- 
sion spectrum and those actually observed will provide 
limits on parameters such as de-hydrogenation, the rel- 
ative abundance of cata-condensed and peri-condensed 
species, and the degree of nitrogen substitution within 
the carb on skeleton, wh ich affects the 6.2 /im C-C stretch 
feature (l Peeters|[2002a|) . 

There is now a great deal of observational evidence 
that PAHs are destroyed within the ionized parts of 
the H II region complexes, with Spitzer observations of 
galactic Hii regions showing clear boundaries between 
the outer PDR P AH emitting zone and the inner pho- 
toionized zone (eg IChurchwell et al.l[200l iPovich et all 
|2007() . The exact destruction mechanism is uncertain, 
but is likely to be photo-destruction through stochas- 
tic heating and/or photoionization and dissociation. To 
simulate this process within the MAPPINGS ill code, 
we previously introduced the Habing photodissociation 
parameter, Ti. — Fpuv/^-HC, a FUV analogy of the stan- 
dard dimensionless ionization parameter U (see eqn 17 
and associated section in SEDl). In a series of test mod- 
els, we found that for typical, solar metallicity starbursts, 
H 10^^ at the ionization front. Henceforth we assume 
this value to be the destruction point for PAHs within 
our models. For typical Hii densities of 10 — 100 cm^'^, 
this im plies a radiation field ~ 2 — 20 times the iHabinel 
(|1968l ) local interstellar radiation field, consistent with 
the range up to which PAHs are observed to survive 
(jCompifegne et al.ll2007l) . At values of H < 10~^, PAHs 
exist in either neutral or singly-charged states, are heated 
by the diffuse FUV/Optical field and emit in the classic 
PAH mid-IR bands. 

This emission spectrum is determined by the natural 
modes of vibration, bending and other deformations of 
the planar carbon skeleton. This spectrum is depen- 
dent upon the size and the electric charge state of the 
molecule, and is modified by the effect of non-hydrogenic 
end groups, including si mple dehydroge nation and skele- 
tal atomic substitution ()Peetersll2002af ). 

Thanks to the advent of space borne IR observato- 
ries, the PAH emission spectrum has now been observed 
in many galaxies and situations , ranging from beautifu l 
maps of Galactic H ii regions (jChurchwell et all l2006f) 



to detail ed mapping of th e features i n both Starburst 
galax ies (|Beirao et al.ll2007[ ) and QSOs (jSchweitzer et al.l 
l2006l) 

Given the wide ranges of possible molecular forms, it is 
surprising that the form of the PAH emission spectrum 
in the Mi d-IR is so siniilar b etween different regions and 
galaxies (|Brandl et al.l [20061) . Only small variations in 
the relative strengths of the PAH features have been ob- 
served within our own Galaxy and in nearby galaxies 
(|Smith et al.ll2007l) . 

The accuracy of using a PAH template to represent the 
series of bands in the mid-IR can be estimated from the 
st udy of the va. r iation of these bands in nearby galaxies 
bv lSmith et all (l2007h . They find on average variations 
of a factor of two around the mean ratios of the differ- 
ent PAH bands (their table 7), with the most signifi- 
cant differences occurring in galaxies hosting weak AGN 
(such as LINERs). This suggests that our PAH template 
is accurate to about this factor, with significant varia- 
tions indicating differences such as PAH ionization state, 
or correspondingly, the presence of a weak AGN in the 
starburst galaxy. The differences between our models 
and the observed PAH bands could therefore be used as 
diagnostics of ISM physics, or host nuclear properties. 

As discussed in SEDl, we parameterize the template 
using a sum of Lorentzian profiles. The Lorentzian fits 
to the spectrum take the form: 



where x — 1/ A cm ^ , the central wavenumber of the 
feature is Xq, the FWHM = 2a, and the peak value is /q 
( ergcm~^s~^Hz~^Sr~^). 

To derive the PAH emission spectrum template cur- 
rently used in these MAPPINGS iii models, we have 
fit Spitzer IRS observations of NGC4676 and NGC7252. 
These two interacting galaxy pairs show strong, clear 
PAH emission, making them good choices for a template 
spectrum. In both objects we subtract the underlying 
dust continuum assuming a combined power law and ex- 
ponential form to fit the PAH-free, long wavelength end 
of the spectra. The combined, continuum-subtracted ob- 
served spectrum is shown with our best fitting template 
in figure [21 The corresponding parameters for each of the 
Lorentzian profiles are given in table [2 

As discussed in SEDl, PAH emission within the MAP- 
PINGS III code is treated as an energy conserving pro- 
cess. In equilibrium all the energy gained by a PAH 
through the absorption of photons can be either lost 
through photoelectric processes or IR emission. Once 
the fraction of the energy lost through photoelectric 
processes is determined using the photoelectric cross- 
sections, the remaining energy fraction is re-emitted in 
the IR according to our empirical template. This is not 
an exact treatment, as the PAH molecules will likely un- 
dergo stochastic heating processes and will lose some of 
their energy through IR continuum emission instead of 
via these fluorescent bands. However the code does allow 
for the stochastic treatment of both very small graphite 
and silicate grains. 

3. THE C PARAMETER 
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We now deal with the parameter which controls pre- 
dominantly the form of the far-IR continuum. Funda- 
mentally, this continuum is a function of the probability 
distribution of the grain temperatures throughout the 
starburst. At any point within an H II region or its sur- 
rounding PDR, this is determined by the intensity and 
the spectral energy distribution of the stellar radiation 
field. Thus, at radius i? in a spherical nebula for any 
individual grain species, s: 

{Tg,.(s)) =/(L,/4^i?2,i?), (8) 

where D is the mean photon energy of the radiation 
field, which depends both upon the age of the cluster 
and the solution of the radiative transfer problem out to 
radius R. Thus, if we wish to find a variable in which 
H II region models evolve along a unique grain temper- 
ature distribution in time, (Tgr(s,f)), these models must 
also preserve the run of L*(t)/47r i?(t)^. Then, because all 
models of this kind give a similar run of grain temper- 
ature distributions, they will also produce very similar 
global far-IR dust emission distributions. With a partic- 
ular choice of cluster luminosity, physically denser H II 
region models have smaller radii, and hence have hotter 
dust temperature distributions. 

We therefore define a Compactness Parameter, C, in 
the form: 

Coci^ (9) 

This variable is akin to the c ompac tness factor found 
in the m odels of iTagaki et al.l (|2003aD and iTagaki et al.l 
(|2003bf) . in that it directly determines the dust grain 
temperature distribution. The ratio on the r ight is also 
comparable to the ratio m,nc/''mc discussed in lSilva et al.l 
(jl998), which controls the SED of their molecular clouds 
and therefore the resulting hot dust SED of their models. 
The stellar luminosity (t) scales with the cluster mass 
Mci, and the radius and pressure at any instant scale 
acc ording to the simp le mass-loss bubble approximation 
of (|Castor et al.lfl975l) : 
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where Lmoch is the instantaneous mechanical luminosity 
of the central stars of the burst (which can be assumed to 
scale as the mass of the cluster) and po the density of the 
ambient medium. Note that the ambient number density 
is no — po/inmn) and ambient pressure Po — nokT^. 
From these equations it follows that. 
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This product remains to be normalized. By analogy 
with the TZ parameter introduced in SED3, we choose 
to adopt the following normalized definition of the Com- 
pactness Parameter; 



logC - log 
5 
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[ Po/k 1 
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+ 5 log 
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(13) 



where Md should now be understood as the mean (lu- 
minosity weighted) cluster mass in the starburst. The 



cm-^Kj ~ 4 and log [Mci/Mq] ~ 3.5 and 
^K] ~ 8 and log [M^i/Mq] - 7, respec- 



likely allowable physical range on log C in starburst envi 
ronments is, roughly, 3 - 7.5. These extremes correspond 
tolog[(P/fc)/ 
log [(P/fc)/cm 
tively. 

In figure[3l we show the variation of with time 

for six different C parameters. These display a strong de- 
crease of the incident flux, and therefore of dust temper- 
ature with time. The specific flux changes several orders 
of magnitude over 10 Myrs. This decrease is due to both 
the increase of the H ii bubble radius with time (see fig. 1 
in SEDl), and the decreasing cluster luminosity as the 
higher mass stars die out. Note that, because of the un- 
derlying power-law behavior of C, on this log scale the 
curves for different C parameters are offset from each 
other in proportion to the change in logC. 

In order to demonstrate the constancy of the dust tem- 
perature distribution at constant logC, in figure 0] we 
show overplotted five SEDs for Solar metallicity pho- 
todissociation regions, which all have the same compact- 
ness parameter of logC = 5.0, but which have different 
pressures and stellar cluster masses. Although these 5 
SEDs appear indistinguishable, they are not exactly the 
same, because their nebular excitation parameters {TZ; 
see SED3) are quite different, and consequently their 
nebular continuum and emission lines are different. 

Provided that we could independently determine TZ 
(from the nebular spectrum) and C (from the form of the 
dust continuum) then in principle we could solve inde- 
pendently for the mean pressure, Po/k and mean cluster 
mass, Mci; 



log 



Mr: 



log 



Po/k - 



'3K 



= logC-l--log7^ 
5 



= log C log TZ. 

5 



(14) 
(15) 



In practice, the separation of these variables would be 
assisted by a direct measurement of the gas pressure. 
For Po/k > 10^ cm~3K we can use the ratio of the [S II] 
A6717, 673lA lines for this purpose. 

To show the effect of varying C on the form of the far-IR 
SED, we present in figure [5] six model PDR SEDs having 
the same metallicity {IZq) and pressure [P/k = 10^), 
with C varied by varying the cluster mass. In the opti- 
cal and near-IR, the model SEDs show stellar emission, 
and the extinction of all six SEDs is the same, as they 
pass through the same column depth of dust and gas. At 
longer wavelengths, the progression in the dust temper- 
atures with increasing C is obvious. 

4. PDR COVERING FRACTION 

In the previous section, our figures |4] and [5] corre- 
sponded to a complete covering fraction of molecular 
clouds; /pDR = 1. In this extreme, the molecular gas and 
dust surrounding the H II regions act as a dust bolome- 
ter, absorbing essentially all of the stellar UV continuum, 
and re-radiating it into the far-IR bump and the PAH fea- 
tures. However, in the case of isolated H II region com- 
plexes in both starburst and in normal disk galaxies, the 
placental molecular cloud is quickly cleared away by the 
stellar winds, and by photo-evaporation. In older clus- 
ters, the disruption of the cluster by this gas ejection will 
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cause the exciting stars to disperse away from the regions 
of high extinction, though the timesc ale for this may 
be gre ater than the H ii region hfetime (jBoilv fc Kroupal 
l2003airb[ ). This process is cutely referred to as "infant 
mortahty" . 

Previously in SEDl and SED2, we parameterized this 
uncovering of the exciting stars by the introduction of a 
molecular cloud clearing or dissipation timescale, Tdoar, 
where the covering fraction of molecular cloud PDRs 
around a stellar cluster is given by; 

/PDR = exp(-t/Tcloar)- (16) 

In SED2, we found that, for Galactic star forming re- 
gions at least, this timescale is quite short, on the order 
of 1 — 2 Myr. However, this certainly does not represent 
all starforming regions, and is probably far too short for 
ULIRGs which have an extremely generous sink of molec- 
ular material. In these objects, the H II regions of indi- 
vidual clusters may merge, but the complex is still sur- 
rounded by molecular gas. Thus, the clearing timescale is 
likely to show a large range, and will depend on the local 
environment. In addition, situations like the commonly 
observed "Blister H ii regions" , where the star formation 
occurs on the edge of a molecular cloud, are not so well 
represented by this formalism. 

To deal with this problem, we introduce here the much 
simpler system of /pdr, which is defined as the time- 
averaged PDR covering fraction during the H II region 
lifetime. Starburst in which the PDRs entirely surround 
their H II regions have /pdr, — 1 while uncovered H II 
region complexes have /pdr — 0. 

In figure [S] we show two series of SEDs to illustrate 
the effect of a changing PDR covering fraction. Both 
represent a solar metallicity starburst with a compact- 
ness C = 10^ and pressure Pa/k — 10^ Kcm~^, and 
are luminosity-weighted integrations scaled to a star for- 
mation rate of 1.0 M0yr~^. The upper set of SEDs 
show the change in SED using /pdr, evolving from a 
pure Hii spectrum (/pdr = 0) to a pure PDR spec- 
trum (/pdr — 1-0). The H II region-only SED (the same 
in both series) is characterized by a strong stellar UV 
continuum, absent or weak PAH features, and hot dust 
emission from within the H II region itself. By contrast, 
the H II region plus PDR models show strongly attenu- 
ated blue and UV continua {Ay ~ 0.8 at /pdr — 1-0), 
strong PAH features and a broad, cool far-IR feature. 

As a comparison, the second plot shows a series in 
7"cicar, goiug from Tcioar = Myrs (pure Hii region) to 
■^cicar = 32 Myrs. Clearly, these two formulations of the 
covering factor are broadly equivalent. On close exami- 
nation it is possible to see some stronger older star fea- 
tures in the Tdear spectra relative to a matching /pdr 
model, but these are small. Concurrent with this is the 
slightly wider IR feature in the /pdr models relative to 
the Tcicar modcls due to the stronger presence of dust 
heated by "old" (~ 10 Myr) star cluster light. 

Note that, as we increase the covering factor, there is 
an increase in the mid-IR from around Afim up to 15/im 
caused by the increasing contribution of PAHs. The far- 
IR dust re-emission feature progressively increases and 
broadens as the contribution of cool dust in the PDR 
becomes more important. The contribution of this cool 
dust in the PDR can be traced at shorter wavelengths 
through the steepening of the 20-35 fim slope. This 



spectral region has few emission lines and there are now 
available ma n y Spi tzer IRS spectra of Starbursts, e.g. 
iBrandl et all |2006[ ). Thus the 20-35 slope may be 
a useful diagnostic of the PDR fraction in starbursts. 
However, this region is also affected by the contribution 
of ultra-compact Hii regions ( §6.1|) and by attenuation 
at high Av {^M- 

Note also the insensitivity of the SED to the covering 
factor in two regions of the spectrum; the 1—4 /im re- 
gion, and at around 20 ^.uy. The 1 — 4 /im emission is 
dominated by the older stars in the starburst, and for 
these the PDR acting by itself is insufficient to produce 
a significant dust attenuation. The constancy of flux 
in the ~ 20/im waveband is somewhat more interesting 
and is directly related to the strong correlation between 
the 24^m Spitzer flux and the SFR (Calzctti ct __al, 2005} . 
The models reveal that almost all of the warm dust emis- 
sion arises from the hot dust embedded within the H ii 
region itself. The cooler dust in the surrounding PDR 
makes little contribution to the global flux at these wave- 
lengths. This result agrees with the spatially-resolved ob- 
servations of nearby galaxies where the 24/im emission 
peaks in Hii regions while the PAH emission is much 
more diffuse. In NGC 5194 about 85% of total galaxy 
24/im emission arises within the defined Hii regions, 
while only ^60% of the 8/im arises within these regions 
(jCalzetti et al.l 120051 ). This result is als o in agreement 
with t he earlier theoretical calculations of lPopescu et al.l 
(|2000[ ) for the edge-on galaxy NGC 891. In this galaxy, 
the star forming disk H II regions have to be associated 
with a dominant hotter dust emission component. 

Our models also demonstrate that the relationship be- 
tween SFR and 24/im emission is not one to one, be- 
cause the warm dust continuum is also sensitive to the 
compactness parameter, C , as is demonstrated in figure 
[5l This finding also agrees with the observations, since 
variations of two or three in the ratio of the 24/im emis- 
sion to the SFR have been observed between galaxies 
(iDale et al.ll200ll: iCalzetti et al.l[2Q0l . 

5. COLUMN DEPTH IN THE PHOTODISSOCIATION 
REGION 

As stated in section 12. 4[ we use a hydrogen column 
density of logiV(HI) — 22.0 (cm~^) to define the extent 
of the photodissociation region. While this value is typi- 
cal of molecular clouds in our own galaxy, it is quite likely 
that those in starburst regions cover a broader range in 
column depths. In figure [7] we show the effect of vary- 
ing the column density in the model. For this exercise, 
we use our fiducial starburst model with solar metallic- 
ity, compactness parameter of C = 10^, ISM pressure 
of Po/fc = 10''^K cm~'^and test a range of PDR column 
depths; logiV(HI) ^ (corresponding to the Hii region 
spectrum), 21.5, 22.0 (our standard PDR model spec- 
trum), 22.5, and 23.0. 

Note that, as more of the UV and visible photons are 
absorbed, the total flux in the far-IR bump becomes cor- 
respondingly greater, and the feature also becomes wider 
as the contribution of the colder dust heated by softer 
photons deep within the PDR becomes more important. 
Note also that the contribution of the PAHs to the SED 
is apparently almost complete by logA^(HI) — 21.5. 

In many respects, the effect of increasing the column 
depth of the PDR is similar to that obtained by applying 
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a diffuse dusty screen (see the following section). As we 
increase in column depth the optical depth increases (rel- 
ative to the H II model) . The effective model Ay for each 
model is 0.2, 0.8, 2.6, and 8.4 as A^(HI) increases from 
to 10^3 cm-2 . At iV(HI) = lO^^cm-^ the IR starts 
to become optically thick, with a corresponding steep- 
ening in the 20/xm to 30/Ltm slope. The increase in the 
silicate feature depth, r9.7^m, becomes obvious somewhat 
earher, at A^(HI) — lO^^cm"^. Note that differences in 
the dust composition and dust geometry ensure that the 
attenuation law of the diffuse dust (shown in fig. [T^ is 
slightly different to that experienced in the PDR, with 
the silicate grains more prominent in the PDR (as seen 
by the stronger silicate absorption feature and flatter ex- 
tinction) . 

The main difference in our model between the PDR 
region and a dusty screen is the inclusion of the self- 
consistent dust emission. By log (HI) — 22 most of 
the UV-optical flux has been absorbed and re-emitted 
as the strong FIR feature. At larger depths the average 
temperature of the dust is quite cold (~ f — 20K) and 
emits only at long (> f00/j,m) wavelengths. However in 
starburst regions such temperatures may not be reached, 
as it is likely that neighboring clusters, as well as the 
underlying diffuse population, would prevent the dust 
from reaching such temperatures. Only in the largest 
molecular clouds, or the most distant dust, would such 
temperatures be reached. It is for this reason that we 
limit our photodissociation regions to logA^(HI) = 22.0 
and use the diffuse population to represent this cool dust 

(El- 

6. OTHER COMPONENTS OF THE SED 

6.1. Ultra- Compact H 11 Regions 

The modeling presented thus far assumes that the 
stars in stellar clusters inflate a common H II region, 
and that the cluster can be treated as a single, cen- 
trally concentrated source of radiation and mechanical 
luminosity. However, in the earliest stages of the clus- 
ter lifetime 10^ years after formation) the individ- 
ual stars composing the clusters are likely to be still 
buried in their separate birth clouds. In addition, we 
know that the massive stars start burning hydrogen even 
before they reach the main sequence and while they are 
still accreting matter from th eir parental molecular cloud 
(|Bernasconi fc Maedeil [TqQGD . During this phase, the 
cluster acts as an ensemble of Ultra- Compact H II regions 
(UCHII) , each trapped at sub-parsec s cales around their 
individual massive parent stars (see IChurchwelll |2002| . 
for a detailed review of UCHII regions). The period of 
time in which the winds of cluster stars cannot operate 
collectively may oc cupy a significant fra ction of a mas- 
sive star's lifetime (jRigbv fc Riekel[2bo4 and references 
within). During this UCHII region phase the optical 
emission will be totally obscured, with ~ 50mag. As 
the UCHI I region is so c ompact, it also displays a hot IR 
emission (jPeetersll2002b( ) . This compactness also ensures 
that the dust is very successf ul in competing aga inst the 
gas for the ionizing photons (|Dopita et al.ll2003l ). which 
makes the region still more compact and ensures that 
radiation pressure effects in the ionized plasma are im- 
portant. Thus, UCHII regions are qualitatively different 
from normal cluster-driven H II regions. 



The technique of modelling the SEDs of individ- 
ual compact and ultra-compact H II regions was 
fully described in 'Dopita^eLaD (|2005bD . Briefly, 
these use the TLUSTY models (available at 
http://tlusty.gsfc.nasa.gov) which were interpolated 
and re-binned to the energy bins used in our code, 
MAPPIN GS III. On the basis of the results presented 
by iMorisset et al.l ()2004D . we can expect the SEDs 
we derive will be very similar to those using either 
the WM-Basic or the CMFGEN models. Unlike the 
TLUSTY atmospheric models, these latter two are fully 
dynamical atmospheric models. The TLUSTY models 
used cover three abundances, 0.5, I.O and 2.QZq. Here 
the definition of solar abundance is effectively the same 
as used in the Starburst99 models. 

We made MAPPINGS iii photo-ionization models for 
the structure of radiation-pressure dominated dusty H II 
regions and their surrounding photo-dissociation regions. 
The models are started close to the central star (1% of 
the computed (dust-free) Stromgren radius) so that the 
ionized gas effectively fills the Stromgren sphere of these 
models. We use the same dust model as for the cluster 
SED computations. The models that we used here have 
fixed external pressure log P/k — 9.0 K cm~'^, and are 
terminated at a H I column density of log A^(HI) = 21.5. 
The mass range of the central stars is 16.7 — 106. 9M0. 
This corresponds to a stellar effective temperature range 
of 32500 < Teff < 52500K, which is covered in bins 
each 2500K wide. The computed SEDs for individual 
UCHII regions correspond to the zero age main sequence 
(ZAMS) of the central stars. 

The starburst galaxy compact H II region SED is con- 
structed by co-adding the SEDs computed for each stellar 
mass bin, luminosity weighted to that corresponding to 
a Kroupa Initial Mass Function (over the effective mass 
range 15 — I20Mq), and scaled to represent a massive 
star formation rate of I.OMq yr~^ continued over a pe- 
riod 0.0- 1.0 Myr. 

The resulting SEDs for each metallicity are shown in 
figure [HI The small changes in the apparent normaliza- 
tion of these spectra is due to the change of the stellar 
luminosity with abundance. Note that all the three spec- 
tra are very similar, with very weak PAH features, a hot 
dust far-IR continuum, similar line spectra and heavy 
(but abundance-dependent) obscuration at shorter wave- 
lengths. 

Because the spectra are quite similar, it is clear that 
any one of them could be used as a UCHII template 
across the full metallicity range. However, for closest 
consistency, we recommend the use of the 2.0.^0 model 
with the 2.0Zq Starburst99 cluster templates, the I.OZq 
model with the I.OZq Starburst99 cluster models and the 
0.5Zq model with the 0.05, 0.2 and 0.4^0 Starburst99 
cluster models. 

As an example, in figure [9] we demonstrate how the 
inclusion of UCHII emission alters the starburst SED. 
We parameterize the addition of this emission to the 
SED by /uchiIj the scale of the UCHII contribution. 
A /ucHii = 1.0 imphes that 50% of the massive stars 
younger than 1.0 Myr are surrounded by ultra-compact 
Hii regions. As the UCHII regions emit predominantly 
in the mid-IR, this parameter affects both the mid-IR 
slope and PAH equivalent widths, as seen in figure [9l 
The emission lines are also affected by the inclusion of 
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the young Hii regions, with the contribution greater in 
the mid-IR than the optical due to the high optical depth 
of the UCHII regions. 

6.2. The Older Stellar Contribution 

The discussion so far has concerned itself only with 
the stars younger than 10 Myr. However, a typi- 
cal starburst will continue for a dynamical timescale 
of a galaxy, typically ~ 10* yr. Over this time pe- 
riod, most of the the young star clusters will have 
dispersed into the field, and away from the active 
star forming regions dSV hitmor e". Chandar fc Falll 120071 : 
iFall. Chandar fc Whitmor&.2005f ). As a consequence, the 
older (> 10 Myrs) stellar population may dominate the 
optical and UV parts of the SED, since it suffers much 
less extinction than the starburst itself, having long es- 
caped its original molecular birth clouds (j Chariot fc Falll 
120001 ). However, it is in the optical-UV range that most 
of the work has been done in constraining the star forma- 
tion history (SFH) of galaxies, using features such as the 
4000A break due to hydrogen and stellar absorption lines 
like the Lick indices to constrain the age and metallicity 
of the dominant stellar population, o r even the evolutio n 
of star formation and metallicity (eg lPanter et al.l[2007t ). 

In order to represent this older stellar contribution in 
our starburst models we use a luminosity weighted sum of 
the 10-100 Myr starburst spectra from Starburst99, hav- 
ing the same parameters as the models used to generate 
the H II spectra. We generate the "old" stellar spectra for 
each metallicity, and assume that the metallicity of the 
starburst the old population are the same. The resulting 
spectra for each metallicity are shown in figure [TUl scaled 
to a continuous star formation rate of l.OM0yr~^. 

The inferred ratio of the starburst fraction in the < 10 
Myr population to the fraction in the age range 10-100 
Myr can provide information about the progress of the 
starburst - whether the starburst activity is accelerat- 
ing or decelerating in recent time. This fraction may be 
constrained by comparing the flux at some wavelength 
in the far-IR bump with the stellar continuum flux at 
any wavelength shorter than about 5/ini. In figure fTTl we 
demonstrate this sensitivity. Here we have added this 
old stellar contribution to a solar metallicity starburst 
with logC = 5 and log P/k = 5. We have normalized the 
SED, assuming a continuous and constant star formation 
rate of 1 MQyr"^ up to the maximum starburst age of 
lO^yr. To emphasize the effect of the older population 
we have taken a totally obscured < 10 Myr starburst 
(/pDR — 1), as shown in the lower curve, and added a 
completely unobscured contribution from the older stars 
to provide the upper curve. 

These models show how it is possible to produce a sys- 
tematic difference between the dust attenuation in the 
UV and the dust attenuation of the H II regions as mea- 
sured by the Balmer Decrement. Such a systematic dif- 
ference has long be known to ex i st from observ ations of 
starbursts (|Calzetti et all 11994 ICalzettil [200l . in the 
sense that the Balmer Decrements indicate greater at- 
tenuation than the UV and visible stellar continuum. 

Finally, the ratio of the two stellar components, 
yo ung/old, is related to the h parameter recently used 
bv iKong et al.l (|2004l ). which is the ratio of the current 
versus past-averaged star formation rate. This parame- 
ter was introduced to help understand the relationship 



between the UV and IR, through concentrating in par- 
ticular on the correl ation of the UV slope and IR excess 
(jMeurer et al.|[l999l ). In our case we parameterize the 
contribution of the old stellar population through the pa- 
rameter /old- Figure [Til presents a case of /om — 1, with 
continuous star formation, and /oid greater and lower 
than one represents cases where the past average star 
formation history is greater and less than the current 
star formation respectively. 

6.3. Diffuse Dust Attenuation 

The simple Hil region plus photodissociation region 
models presented in the previous sections provide a max- 
imum attenuation of Ay ~ 0.8 for a solar metallicity 
starburst and Ay ~ 2.5 fo r 2Z[7^. This i s low er than ob- 
served in many ULIRGs (|Farrahet al.ll200l . In order 
to properly account for the total attenuation, we need to 
include the attenuation produced by a by a foreground 
dusty screen associated with gas in the starburst host 
galaxy, but not necessarily partaking in the starburst, 
to account for these heavily obscured starbursts. This 
foreground screen will also attenuate the older stars as- 
sociated with the starburst (discussed in the previous 
section). 

The properties of such a turbulent foreground at- 
tenuating screen were discussed in a series of papers 
by F ischera & Dopita (jFischera fc Dopital 1200^ 12004 
|2005[) . They showed that turbulence which produces a 
log:normal distribution in local density will also, to a high 
level of approximation, produce a log:normal distribution 
in column density. They also showed that the resulting 
attenuation curve is unlike that of a normal extinction 
law, showing lower attenuation in the UV and larger at- 
tenuation in the IR, due to the spatially-varying extinc- 
tion across the face of the dust-obscured object. Here 
we adopt the theoretical attenuation curve computed in 
those papers, which is shown in fig. [12] and which pro- 
vides a close approximation to the empiricall y-derived 
Calze tti extinction law for starburst galaxies (|Calzettil 
[200l . This curve does not allow for the possible destruc- 
tion of PAHs in the diffuse medium, and the computed 
2175 A absorption feature may be rather too strong to 
be applicable to starburst environments. 

In figure[T3|we show the effects of this dusty screen on a 
starburst with solar metallicity, compactness parameter 
of C = 10^ ISM pressure of Po/fc = lO^K cm'^, and 
covering fraction of /pdr = 0.5. We show both low and 
high attenuation, with 10 SEDs plotted in figure [T3] with 
Av of 0.0, 0.5, 1.0, 2.0, 4.0, 6.0, 10.0, 15.0, 30.0 and 
40.0, with a clear depletion of the ultraviolet and optical 
flux as the Ay increases. To emphasize the effects of 
the attenuation on the IR emission we do not include 
the diffuse cool-dust emission that would be expected 
with such attenuation. We note that beyond ~ 60/im 
there would be a contribution due to thermal emission 
from this cool dust, and care should be taken with any 
interpretations made beyond this wavelength. The effect 
of this emission is discussed in the following subsection. 

The attenuation and reddening of the underlying 
starlight at the optical and UV wavelengths is very evi- 
dent. However, it is only at Ay <; 5 that absorption in 
the 9.7 /im and 18 /im silicate featu res becomes appar- 
ent. With Av/T9.7^m ~ 16.6 (,Rieke fc Lebofskvl 119851) 
this feature is weaker than the optical opacity, but it 
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is observed to be visible or very strong in a large num- 
ber of starburst galaxies, proving that a good deal of 
the starburst activity will be totally obscured in the vis- 
ible thanks to large column densities in the surrounding 
molecular gas. Note also that at the highest Ay even the 
20-35/zm slope steepens to the "reddening" effect of the 
dust attenuation. 

The Ay > 20mag. dusty screen needed to provide the 
observed depth of the silicate absorption troughs seen 
in some ULIRGs has little effect on the global energy 
balance of the starburst. By an Ay of 1.0, most of the 
FUV-optical radiation from stars, almost 80%, has been 
absorbed (or scattered) by dust. The dechne in the lumi- 
nous flux through the PDR, accompanied with the soft- 
ening of the radiation field naturally produc e s the dust 
temperature gradient which iLevenson et al.l (|2007| ) be- 
lieve is required in the obscuring material. 

6.4. Diffuse Dust Emission & Scattering 

A complete model of a starburst environment should 
not neglect the contribution to the far-IR emission by 
diffuse cool dust. Some of this may well be the same dust 
that produces the foreground screen absorption. The 
optical photons that are absorbed by the diffuse galactic 
dust at high Ay are still capable of heating the dust in the 
diffuse ISM. In addition, the star forming regions may not 
be fully cocooned by its surrounding PDR cloud, so that 
some fraction of the cluster UV radiation may escape and 
heat this diffuse dust. This cool galactic dust will mostly 
contribute to the SED in the > 100/im to sub-mm region. 
A portion of this radiation may also be scattered into our 
sight-line, which mostly affects the far-UV SED. In our 
modeling we have not included this component due to 
our limited geometry. 

The exact temperature distribution of the grains will 
also depend upon the geometry of the starburst, pre- 
existing stellar population, and dust within the star- 
bursting galaxy. Due to the limited geometry of our 
simple models, we are unable to m odel such d i stribu - 
tions. Instead we follow the work of iDale et all (|200lh . 
and calculate the diffuse cool dust emission in terms of 
the average interstellar radiation density. This allows us 
to model both distributed starbursts where the average 
radiation field is high, and nuclear starbursts, where the 
rest of the galactic scale dust is only weakly heated. 

To represent the average radiation field we have used 
the luminosity-weighted average of a Starburst99 cluster 
from 10 — 100 Myr discussed in the section 16.21 While 
this is younger than the average age of the pre-existing 
stars in a starburst host galaxy, these stars are nonethe- 
less likely to provide the dominant dust-heating radiation 
field. 

We then calculate the cool dust emission using the 
MAPPINGS III code, assuming the same dust proper- 
ties as before and a column depth of log N{H) = 22.5 
( cm~^). For the he ating flux we scale the radiation field 
in terms of the local iHabind (j 19681 interstellar radiation 
field (ISRF: FUV~ 1.6 x lO'^ergs'^ cm'^). We set the 
radiation field from 0.1 to 1000.0 times this value, in steps 
of 0.3 dex. The resulting IR emission is shown in figure 
fT4l As expected, low ISRF leads to very cold dust, and 
the high ISRF has dust temperatures similar to those 
encountered in our PDRs. 

The actual contribution of this diffuse dust emission 



component to the global starburst SED is not con- 
strained by these models. This would require a more 
sophisticated geometrical model of the starburst, its out- 
flows, and any more extended disk or tidal structure 
around the starburst core. 

In Figure [15] we show the effect of adding this diffuse 
dust emission component with a total intensity of 10% of 
our fiducial model Starburst with solar metallicity. For 
clarity we add only four of the diffuse emission mod- 
els, heated by a Habing radiation field, Gq of 1000.0, 
100.0, 10.0, and 1.0 times the local interstellar radiation 
field (ISRFiocai), respectively. The attenuation associ- 
ated with the diffuse dust emission is not included here. 

The model with 1000.0 times ISRFiocai has hotter dust 
than the log C=5 PDR and is therefore likely to be un- 
physical. The 100.0 times ISRFiocai model has a sim- 
ilar dust temperature as the PDRs, and so the diffuse 
field serves only to increase the total flux of the IR. This 
probably represents the extreme case for a starbursting 
galaxy, and may lead to the narrow and strong far-IR 
features seen in some ULIRGs. 

The cooler diffuse emission models, with Go < 10.0 
times ISRFiocai, are both applicable to less energetic 
starburst galaxies. In such cases, the diffuse emission 
acts to broaden the IR feature, as well as shift the peak 
to longer wavelengths, while leaving the shorter wave- 
lengths (< 60/im) relatively unaffected. 

7. STARBURST METALLICITY 

The metallicity of a starburst affects the SED in sev- 
eral ways; through the intrinsic change in the stellar SED 
with metallicity, through the changing gas-phase abun- 
dances, which determine the temperature and the line 
emission of the H II regions, through the opacity of the 
ISM in the dust, and through the metallicity-dependent 
change in grain composition. A full discussion of the 
effect of metallicity on the emission line spectra of the 
H II regions was given in SED3. In this section we will 
systematically investigate the remaining effects. 

In figure [16] we show our fiducial (H II region & PDR) 
starburst with log C = 5 and \ogPo/k = 5 computed 
using the 5 standard Starburst99 metallicities. As the 
PDR is defined through a constant column depth of hy- 
drogen, lower metallicity leads to lower column of dust. 
This leads to a strong decrease in the optical-ultraviolet 
opacity as the metallicity is decreased. 

The metallicity-dependent effects on the H II emission 
line spectrum have previously been remarked upon in 
SED3. As the metalHcity decreases, the stellar spectrum 
becomes harder due to the decreasing opacity of the stel- 
lar atmospheres and winds. In addition, for a given size, 
a stellar cluster of lower metallicity has a higher ionizing 
luminosity and lower mechanical luminosity. This leads 
to more compact H II regions, and higher ionization pa- 
rameters in the surrounding H II region. The fraction of 
radiation absorbed by the dust in the H II region de- 
pends o n the product of m etallicity and ionization pa- 
rameter (jPopita et al.l[2003h . This product remains ap- 
proximately constant with metallicity, ensuring that the 
flux under the far-IR bump remains approximately con- 
stant for the H II region SEDs. 

In the PDRs, the UV to optical SEDs clearly shows the 
effect of increasing opacity with metallicity. The far-IR 
features change in several ways. First, there is a system- 
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atic increase in the relative strength of the PAH emission 
features with metalhcity. This is a consequence of the in- 
crease in the C/0 ratio with metalhcity, which ensures 
that PAHs account for more absorption at higher metal- 
hcity relative to the silicates, combined with the higher 
mean dust temperatures which characterize lower metal- 
licities. Another factor is the increased strength and 
hardness of the average radiation field in the constant col- 
umn PDR with decreased metalhcity. As a consequence, 
the Habing PAH survival criterion is only met deeper in 
the cloud, resulting in overall weaker PAH emission. 

Such a decrease in the PAH strength with decreas- 
ing metalhcity has been observed with both the Spitzer 
and ISO Space Observatories jE ngclbracht ct al. 200^ 
Rosenbers et al. 2006; JWu et all I2006t iMadden et ail 
2006 : Jackson c t al. 200^ However, the observed de- 
pletion of PAHs in the low-metallicity environments may 
be ev en greater than that c omputed in ou r mod els, and 
both iWu et all (|2006f ) and iJackson et al] (|2007D impli- 
cate grain destruction processes as acting more efficiently 
in the lower metalhcity environments. We would hope 
that our models, applied to these low-metallicity systems, 
would be able to better confirm and quantify this effect. 

In the PDRs, the IR flux can also be seen to increase 
and become broader with metalhcity up till Z ~ I.OZq, 
after which it stays approximately constant. This is pre- 
dominantly due to the increased dust column. The shift 
of the IR peak to longer wavelengths is also partly due to 
the increased dust column but is also due to the increas- 
ing mechanical luminosity of the starburst with metal- 
hcity, which results in larger Hii regions, with cooler 
average dust temperatures. 

8. COMPARISON WITH OBSERVATIONS 

8.1. Data Sources 

In order to compare the models with data, we require 
as close a homogeneous data-set as available, covering as 
wide a wavelength range as possible. For this purpose, 
we have selected a pair of popular template starburst 
galaxies, NGC 6240 a nd Arp 220 from the 41 ULIRGS 
observed with ISO bv lKlaas et~al] JlOOl). This data set 
is ideal for our purpose because their SEDs are well sam- 
pled over the full wavelength regime 1-200 /im. 

Flux densities at other wavelengths were collected us- 
ing the NASA/IPAC Extragalactic Database (NED) sup- 
plemented with a wide selection of on-line catalogues 
and papers. Generally, the UV/optical fluxes are taken 
from the third reference cat alogue of bright galaxies v3.9 
(jdeVaucouleurs et al.lll99l|). Manv opt i cal an d near-IR 
fluxes were taken from ISpinoglio et al.l (|1995f ). or from 
the APM and 2MASS databases. 

The majority of data poi nts in the 1- 1 300 ju m wave - 
length range come from iKlaas et all (Il997l. |2001D 



Additional data points w ere taken from iSanders et al 
(•igsM), 'Sander s et all (ll988bD . Murph v et al 
1996|) . iRigopoul ou et al.l (I1996D. [Rig opoulo u et al 
1993) . ISu^ace fcSandersI (|2000|). [Lisenfeld et al. 
2000), lDunneetal.1 (|2000f ) Jpunne fc EalesI (|2001h . 



IScoville et all (|2000f ). and ISpoon et all ([2004) and 



to near-IR fluxes have been corrected for Galactic extinc- 
tion using the E(B— V) values based upon IRAS 100 /xm 
cirrus emission maps ([Sc hlegel e t al.lll998l ) and extrapo- 
lating following ICardelir et al [8^. 

For comparison of the models with a typical Spitzer 
IRS spectrum of a starburst, we have used the latest (re- 
calibrated) v e rsion of the spectrum of NGC 7714 from 
iBrandl et all (|2006h . 

8.2. Pan-Spectral Fitting 
Our library of models contains the following elements: 

• An ensemble of H II regions surrounding young 
clusters with ages < 10 Myr characterized by mean 
compactness C and metalhcity Z . 

• A set of PDRs surrounding these H II regions with 
a mean geometrical covering factor, /pdr - 

• A population of young (< 1.0 Myr) ultra-compact 
H II regions and their PDRs surrounding individual 
massive stars, characterized by a fraction /uchii, 
where /uchii = 1.0 would imply that 50% of mas- 
sive stars younger than 1.0 Myr were surrounded 
by ultra-compact H II regions. 

• An older stellar population with ages 10 < t < 100 
Myr. The flux of this component is scaled by a 
factor /old, where /om = 1 would correspond to 
continuous star formation over a total period of 
100 Myr. 

• A foreground turbulent attenuating dust screen, 
characterized by an optical depth in the V- band. 

Ay. 

• A re-emission component from the diffuse ISM, 
characterized by the mean Habing fleld intensity 
Go and scaled to a percentage of the total Bolo- 
metric Flux of the Starburst (< 20%). 

Ideally, all of these elements should be fltted via a non- 
linear least-squares procedure. However, lacking this tool 
for the time being, we have elected to make hand-crafted 
flts to a few template objects. In addition, we have cho- 
sen to make the following simplifying assumptions: 

• We treat the < 10 Myr stellar population as an 
ensemble of H II -|-PDR regions with /pdr = 1-0. 
This approximation is justifled by the observation 
that the global obscuration of a starbur st increases 
with t h e absolute rate of star f ormat i on iBuat et al 



19991) : lAdelberger fc St eidel (|2000l ): iDopita et al 
20021 ): IViih et all (20031. This is a natural conse- 



references therein. 

When available, the UV/optical and near-IR {JHK- 
band) points include aperture corrections to allow a di- 
rect comparison with th e larger aperture mid- and far-IR 
fluxes (see, for example. ISpinoglio et al.l (jl995l )). All UV 



quence of the lKennicuttI ((19^) star formation law, 
connecting the surface density of star formation 
to the surface density of gas; SgpR cx: Sg^s, with 
n ~ 1.3 — 1.6. Noting that Ay oc Sgas, it follows 
that intense starbursts are highly dust-enshrouded. 

• We ignore the contribution of the diffuse dust in 
the far-IR emission. As has already been noted in 
Section 16. 3[ the contribution this makes is small, 
and significant only above 100 /im. 
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With these assumptions, we show in figure [17] the fit 
we obtain for the galaxy NGC 6240. As can be seen, 
with five free parameters we obtain an excellent fit to 
the observations over nearly 3.5 decades of frequency. 
In this figure, the data have been scaled to a star for- 
mation rate of 1.0 M0yr~^ in order to make this SED 
comparable to the others in this paper. The scaling fac- 
tor required indicates a total star formation rate of 120 
MQyr~^ for this galaxy, a little higher than the value 
of 102 M(7iyr~^ derived from the IRAS far-IR luminosity 
bv lDopita et al.l (|2002f ). However, the Ha luminosity f or 
this galaxy (also measured by bv iDopita et al.l (|2002t l). 
which comes from the extended gas indicates a star for- 
mation of 14.6 MQyT~'^ , provided that it is dominated 
by star formation, and not by an obscured active nucleus 
(see below). Since the far-IR comes from the obscured 
starburst, and the visible Ha from the regions which are 
relatively unobscured, these figures suggest that the total 
star formation rate is indeed close to 120 M0yr~^. 

As can be seen in figure fTTl with a model with only five 
free parameters we have obtained an excellent fit to the 
observations over nearly 3.5 decades of frequency. How- 
ever, this quality of fit to a model which only includes 
elements of a starburst system is at first sight extraor- 
dinary. NGC 6240 has long been implicated with an ac- 
tive nucleus, and shows both a str ong non-thermal radio 
exces s, extended radio emission (jCallimore fc Beswickl 
l2004f ). and X-ray e mission as sociated with a highly dust- 
obscured AGN (egJ IkebeeLal. 2000; Risaliti et aij |2000l : 
iKewlev et all 120001 ) . Recently lArmus et al.l (|2006l ) have 
directly detected the active nucleus via the [Ne V] 14.3 
/zm line using the IRS on the Spitzer Space Telescope. 
From this, they estimate that the AGN has a flux of 
3 — 5% of the bolometric luminosity. At such low levels, 
it would account for the slight excess in the fiux in the 
vicinity of the 6 — 14//m PAH features seen in the obser- 
vations when compared to our model, but is not sufficient 
to compromise the rest of the fit. 

NGC 6240 has an unusually strong contribution from 
older stars, suggesting either that the starburst is rather 
old in this object, or that the current starburst is the 
second episode in this galaxy, or that there is an impor- 
tant contribution of stars older than ~ lO^yrs. All of 
these hypotheses are consistent with the known status of 
NGC 6240 as a post-merger system. 

What about the uniqueness of this fit? Fortunately, 
the different parameters of the fit act on different parts 
of the pan-spectral SED, so it is fairly easy to separate 
them when we have data covering such a large range in 
wavelength. We have performed the test of varying each 
of the major parameters, excluding metallicity, system- 
atically around the best-fit solution, and the result is 
shown in figure [TSl 

The effect of varying metallicity is shown in figure fT6| 
and has perhaps the most profound effect on the spec- 
trum. To summarize, the metallicity is most easily de- 
termined from the shape and absorption line intensities 
of the stellar continuum, and by the emission line tech- 
niques discussed in SED3. It also has an effect on both 
the strength of the PAH features and the width of the far- 
IR bump. High metallicity gives strong PAH features and 
wide far-R bump. In the fitting described in this section 
we have mostly relied upon these latter characteristics to 
estimate the abundance. 



For the remaining parameters, figure [18] shows that 
each affects a different part of the SED. A change in logC 
only changes the far-IR bump, shifting it in peak wave- 
length without appreciable change in the total width. 
A change in /om simply scales the 0.091 — 5/im spec- 
trum up and down, leaving the rest of the SED un- 
changed. A change in Ay affects the slope of the visible- 
UV spectrum. Note that, when Ay is higher than about 
5 — 10 mag, the 9.7 and 18.0/im silicate absorption fea- 
tures appear and can be used to constrain the extinc- 
tion when the the visible-UV part of the spectrum is 
too attenuated. The Ultra-Compact H II region fraction 
/ucHii mostly affects the 10 — 30^m part of the spec- 
trum. In figure [THJ the apparent sensitivity of the SED 
to this parameter seems small, but this is because the 
inferred logC for this galaxy is very high. Thus, the fiux 
of the ordinary H II regions is high in the 10 — lOO/im 
region where the UCHII population is important, and 
the contribution of the UCHII regions is veiled. Galax- 
ies with lower logC show the UCHII contribution much 
more clearly (see fig. [S]). 

As another example of a famous starburst, we present 
in figure [19] our fit to Arp 220, the predominantly used 
template of starburst SEDs. This galaxy is character- 
ized by a greater optical extinction, lower compactness 
parameter, and a lower contribution of the old star pop- 
ulation than NGC 6240. In addition, a lower chemical 
abundance is indicated by the weaker PAH features, and 
the narrower, more sharply peaked far-IR bump. 

Arp 220, as the local ULIRG, has quite often been used 
as a testb ed for SED mode lling. Examples of these can 
be seen in ISilva et all (|l998l the ir figure 9).lTagaki et all 
l 2003bL their fi gure 8) and ISiebenmorgen fc Kriigell 
1 2007L their figure 5) to name a few. While all these 
models (including our own) can be seen to fit quite rea- 
sonably the available observations of Arp 220, they differ 
in their model parameters, making direct comparisons 
difficult. However the main physical conclusions drawn 
from the models are the same, and comparisons here can 
give insight into both the models and Arp 220 itself. 

All models suggest a star formation rate (SFR) for 
Arp 220 of '^300 M0yr~^. We derive a star formation 
rate of 315 M0yr~^, w hich can be c ompare d with the 
270 M (nyr~-^ obtained by Shiova. Tren tham & Taniguchi] 
(|2001h . 260 M^yr-i by iTagaki et all (t2003b|) and 580 
Mgyr-i bv ISilva et all mm . 

In connection with this is the total luminosity of 
Arp 220, which is logL* = 12.16(Lm) in our mo dels, 
close to the value of 12 . 1 of iTagaki et al.l (|2003bD and 
ISieben morgen fc Kriigell ()2007| ). and just below the 12.4 
from [Silva et al. (1998). For the total stellar mass of 
Arp 220 we obtain logM* ~10.5 (M0), higher than the 
previous SED model estimat es of 10.4 (,Silva et al.lll998D 
and 10 (jTagaki et al.ll2003bh . but due to the low mass- 
to-light ratio of our older stellar component ( §6.2p this 
value is somewhat uncertain. All these values indicate 
that both our model and fit to Arp 220 are at least con- 
sistent with previous models, and suggest the true values 
for Arp 220. 

One well known detail of Arp 220 is its 

very high nuclear extinct ion; Ay ^ 30 

(IShiova. Trentham fc Taniguchil 120011 : ISpoon et all 
l2004f ) has been e stimated and even higher est imat es 
from models exist (jSiebenmorgen fc Kriigell [2007| ). This 
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illustrates a limitation of our fitting procedure. The 
derived Ay is determined essentially from fitting the 
attenuation of the older stellar population, not the 
nuclear region. When both the 9.7 and 18.0/im silicate 
absorption features are observed along with optical or 
near-IR stellar continuum, we could then modify our 
fitting procedure to first fit the Ay of a foreground 
screen for the H II regions + PDRs implied by the depth 
of the silicate absorption, and then apply a second, more 
optically-thin foreground screen to match the extinction 
of the older stars seen in the visible-UV part of the 
spectrum. 

In order to test the effect on the fitting parameters 
we have made such a model, which we present in figure 
[20l This model, with Ay = 20 mag for the starburst 
component and Ay = 2.2 for the older stars certainly 
fits the region of the silicate absorption features better, 
and allows a much larger fraction of compact H II regions; 
/uCHii 0.7 - implying that ~ 41% of all stars younger 
than 1.0 Myr are found in compact H II regions. 

The fact that we derive a lower Ay for the starburst 
than other authors is not surprising, since the attenua- 
tion law that we use provides greater attenuation in the 
IR and less in the UV than a standard extinction law, 
thanks to the patchy nature of the foreground screen. 
What is surprising is the reduction in the strength of 
the emission lines in the visible and UV regions of the 
spectrum. This is caused by the much larger nuclear at- 
tenuation, and shows that the equivalent widths of the 
IR emission lines compared with the visible or near IR 
emission lines may be used as a sensitive diagnostic of 
nuclear extinction. 

It should be noted that, while Arp 220 is one of 
the predominantly used starburst templates, there is 
increasing evidence that this object is not a good 
representative for high-z starfor ming galaxies (see eg. 
iMenendez-Delmestre et al]l2007f ). Rather, less extreme 
objects such as M82 or NGC 6240 are better local analogs 
of the high-z actively star-forming objects such as sub- 
millimeter galaxies. 

8.3. Fitting Spitzer IRS Spectra 

As a final example of the fitting process, we compare 
our fit with the detailed Spitzer Space Obs ervatory IRS 
low resolution spectra of NGC 7714 from iBrandl et al.l 
((2006). The fit is shown in figure [HI Note that this ob- 
ject has a much lower star formation rate (8.0 Moyr"^) 
than the previous two examples. The star formation rate 
is very well constrained by the normalization process, 
and can be determined to an accuracy of better than 
5%, assuming that the IRS aperture integrates the full 
extent of the star forming region. 

For this object, the older stellar component is not well- 
constrained, as the IRS spectra do not quite extend to 
short enough wavelengths to measure it. Also, because 
the attenuation is not enough to produce appreciable 
9.7/im silicate absorption {Ay < 5 mag.), we cannot mea- 
sure Ay from these spectra, so it has been set equal to 
zero. Note that the line emission spectrum is quite well 
fitted by the model, except for the strength of both the 
[S IV] and the [Ar III] lines in the vicinity of the 9.7/im sil- 
icate absorption band. Again, this may indicate a rather 
higher obscuration of the young H II regions than in the 
model. 



In this spectrum the abundance is fairly well con- 
strained by the strength of the PAH features and the 
equivalent width of the emission lines, while the 20 — 
35/im slope constrains the values of logC and /uchii- 
The rather steep slope in the continuum spectrum at 
about 16/im is a characteristic signature of the presence 
of compact H II regions. 

9. DISCUSSION & CONCLUSIONS 

In this paper we have described an extensive library of 
pan-spectral SED models applicable to starburst galax- 
ies, and demonstrated the promise of these in deriving 
the physical parameters of starbursts. These models rely 
upon a local, rather than a global solution to the radia- 
tive transfer. Such an approach works because of the fact 
that, in starburst galaxies, the vast majority of the far- 
IR emission arises from absorption of the UV radiation 
field in a relatively thin dust layer, the classical photo- 
dissociation region (PDR). This region has a typical op- 
tical depth corresponding to Ay ~ 3, and a thickness 
AR ^ 300/nH pc. In the molecular regions surrounding 
normal galactic H II regions, hydrogen densities are typ- 
ically 100 — 1000 cm~'^, implying that much of the far-IR 
they produce comes from a layer of parsec or sub-parsec 
dimensions. In starburst galaxies, interstellar pressures 
may range up to a factor of 100 higher than this, pro- 
ducing correspondingly thinner PDR zones. In addition, 
the Stromgren volume, the volume of the ionized gas in 
the H II region surrounding the exciting star or cluster 
scales as n^^, making H II regions much more compact 
as the pressure in the ISM is increased. 

By simplifying the radiative transfer problem to a lo- 
cal one connected with individual clusters and their H II 
regions and PDRs, we can compute the SED as the sum 
of a set of effectively independent components. Our li- 
brary of models provides the following ingredients to the 
pan-spectral SED of starbursts: 

• An ensemble of H II regions surrounding young 
clusters with ages < 10 Myr. 

• A set of PDRs surrounding these H II regions. 

• A population of young (< 1.0 Myr) ultra-compact 
H II regions and their PDRs surrounding individual 
massive stars. 

• An older stellar population with ages 10 < t < 100 
Myr. 

• A foreground turbulent attenuating dust screen. 
Separate screens may be used for the younger < 10 
Myr population and the older stellar population. 

• A re-emission component from the diffuse ISM. 

We have shown that the position of the far-IR dust 
re-emission peak is primarily controlled by Compactness 
Parameter C defined in section [3l although the position 
and shape of this feature is also influenced by the mean 
covering fraction of the PDRs surrounding the individual 
H II regions, /pdr, investigated in Section 01 and by 
the metallicity discussed in Section [T) In addition we 
have investigated the effect of the column density in the 
PDRs surrounding the H II regions, provided a global 
spectrum of an ensemble of compact H II regions derived 
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from our earlier work, and have investigated the effect 
of metalhcity on the SED of the older stars with ages 
10<t< 100 Myr. 

Finally, we have provided the attenuation properties of 
a turbulent absorbing dusty screen, and have computed 
simple one-dimensional models of the thermal emission 
from the diffuse dust illuminated by the SED of the older 
(> 10 Myr) population. These models are characterized 
by the local diffuse radiation field intensity, expressed in 
units of the Habing field intensity. Go, where Gq = 1.0 
corresponds to the intensity of the diffuse radiation field 
in the vicinity of the sun. 

We have demonstrated how the far-IR to sub- mm SED 
is controlled by the Compactness parameter logC, and by 
the metallicity. This is of particular applic ation to the 
high-r edshift sub- mm galaxies. As shown by lBlain et al.l 
()2004( ). a modified Black-Body fit can be made to the 
long wavelength side of the far-IR peak in starburst 
galaxies to derive a "dust temperature". Although the 
concept of such a dust temperature is physically mean- 
ingless in the light of our models, which contain a wide 
distribution in dust temperatures, it is a useful way to 
characterize the slope and the position of the sub-mm 
SED, and may well be related to the minimum dust tem- 
pera ture in th e starburst. 

Bl ain et all £2004) showed that, in ULIRGs, the "dust 
temperature" derived in this way is observed to corre- 
late with the absolute luminosity (or equivalently, to the 
rate of star formation). In our interpretation we would 
conclude that for ULIRGs, the compactness parameter 
increases with increasing luminosity. This would be con- 
sistent with more luminous galaxies having greater sur- 
face densities of star formation and greater gas pressures 
and densit i es. Th is in turn is in accord with the empirical 
iKennicuttI (|199 8D law of star formation; Ssfr oc Sgat=^° ^ 

However, the lBlain et al.l (2004 ) work also showed that 
the high-redshift submillimeter selected galaxies (SMGs) 
provide a similar correlation, but shifted to higher lu- 
minosity. At a given luminosity, the dust "tempera- 
ture" in SMGs is about 20K cooler than in ULIRGs in 
the local universe, and at a given dust temperature, the 
SMGs are typically 30 times as luminous as their ULIRG 
counterparts. Given that we have no reason to suspect 
lower dust temperatures in sub-mm galaxies, we must 
conclude that the starbursts in these galaxies have simi- 
lar compactness to local starbursts, but are typically 30 
times more luminous and spatially-extended than local 
ULI RGs. 

Ta gaki et al.l (|2003al lbl) had previously found that most 
ULIRGS have a constant surface brightness of order 
lO^^LQkpc"^. These parameters probably characterize 
"maximal" star formation, above which gas is blown out 
into the halo of the galaxy and star formation quenched. 
In order to scale the star formation up to the rates in- 
ferred for SMGs (~ 1000-5000Moyr-i), we must involve 
a greater area of the galaxy in star formation, rather 
than trying to cram more star formation into the same 
volume. For a typical value of lO^^L0kpc^^, we require 
"maximal" star formation over an area of ~ lOkpc^, and 
the most luminous SMGs require star formation to be 
extended over an area of at order ~ lOOkpc^. 

From our library of models we have constructed syn- 
thetic pan-spectral SEDs composed of only the following 



components: 

• The young H II regions with ages < 10 Myr charac- 
terized by a mean compactness, C and metallicity. 

• The PDRs surrounding these H II regions with a 
mean geometrical covering factor, /pdr- 

• A population of young (< 1.0 Myr) ultra-compact 
H II regions and their PDRs surrounding individual 
massive stars, characterized by a fraction /uchii, 
where /uchii = 1-0 would imply that 50% of mas- 
sive stars younger than 1.0 Myr were surrounded 
by ultra-compact H II regions. 

• The older stellar population with ages 10 < t < 
100 Myr, scaled by a factor /om, where /om = 1 
corresponds to continuous star formation over the 
period of 10 < i < 100 Myr. 

• A one- or two- component foreground turbulent at- 
tenuating dust screen. 

Together, these components allow us to construct theo- 
retical pan-spectral SEDs which encompass the full rich- 
ness and variety observed in the SEDs of real starburst 
galaxies. We note that different parts of the SED are 
sensitive in different ways to these various theoretical 
components. Therefore, given a sufficient spectral range 
in the observations, we can create fits to observed SEDs 
which are unique, and which allow us to extract each 
of the various physical parameters listed above for the 
individual starburst galaxies. Finally, the scaling of the 
bolometric flux of the theoretical SED to the observed 
bolometric flux of the starburst provides us with an ac- 
curate estimate of the total star formation rate , mea- 
sured in M0yr~-'^. 

We have demonstrated the utility of this method by 
fitting theoretical SEDs to a few famous template star- 
bursts such as Arp 220 and NGC 6240. These fits have 
been "hand-crafted" , and it remains to automate the 
process to a multi-variate least-squares fitting procedure, 
which we hope to present in a future paper. 

Although this paper has dealt exclusively with star- 
bursts, we note that many Ultra-Luminous Infrared 
Galaxies (ULIRGs) are known to contain an active galac- 
tic nucleus. To deal with these cases, we need to add a 
further five components to the mix: 

• The UV-IR continuum emission of the AGN itself. 

• Hot dust from the accretion disk around the AGN, 
emitting mostly in the 5 — 50/im waveband. 

• The global emission from the extended Narrow 
Line Region (NLR) surrounding the AGN. 

• A foreground dust screen surrounding the accretion 
disk, the so called "dusty torus" 

• Non-thermal radio synchrotron emission compo- 
nent or components. 

This will be the subject of a future paper in this se- 
ries, but we note that many of the required compo- 
nents are already fairly well understood. The AGN 
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itself is usually approximated by a power-law, or bro- 
ken power law, a nd the NLR c omponent has already 
been computed by lGrove"i (|2006D . The radiative trans- 
fer through a dusty torus around the nucleus is a 
much more complicated issue that has been dealt with 
by several authors. Originally pictured and mod- 
elled as a smooth dis t ribution (iPier fc KrolikI [l99^ 
Granato &: Danes^ll994 lEfstathiou fc Rowan-RobinsonI 



1995D . more recent models assume a more clumpy 



structure to explain the observed distribution of the 
10/im silicate f eatur e and width of the far-IR peak 
(iNenkova et all 120021 : iDullemond fc van Bemmel 120051 : 
iHonig. et al.ll200l . In any case, treatments of the vari- 
ous components of an AGN SED do exist, and must be 
considered when fitting the SEDs of ULIRGs and high- 
redshift luminous objects. 

The replacement the of semi-empirical approach by the 
new quantitative approach to SED fitting of starbursts 
presented here should greatly enhance the utility of exist- 
ing and future surveys, by providing detailed estimates 
of the physical parameters of the starbursts. In this way, 
it should assist in the derivation of statistical parame- 
ters, demography and cosmic evolution of this class of 



object, and should cast light on the nature of starbursts 
in the high redshift universe, including sub-mm galaxies 
and the high-redshift radio galaxies, when today's mas- 
sive "red-and-dead" Elliptical galaxies were in the pro- 
cess of assembly, and ULIRGs ruled the star formation 
rate density of the Universe. 
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TABLE 1 
The Solar Abundance Set 
(Zq) and logarithmic 

DEPLETION FACTORS LOG(D) 
adopted for EACH ELEMENT. 



Element log(Z0) log(D) 



H 


0.00 


0.00 


He 


-1.01 


0.00 


C 


-3.59 


-0.52 


N 


-4.20 


-0.22 


O 


-3.34 


-0.22 


Ne 


-3.91 


0.00 


Na 


-5.75 


-0.60 


Mg 


-4.42 


-0.70 


Al 


-5.61 


-1.70 


Si 


-4.49 


-1.00 


S 


-4.79 


-0.22 


CI 


-6.40 


-0.30 


Ar 


-5.44 


0.00 


Ca 


-5.64 


-2.52 


Fe 


-4.55 


-2.00 


Ni 


-5.68 


-1.40 



TABLE 2 
Normalized Parameters 

OF the Lorentzian 
Components of the PAH 
Emission Band 











al 


3040.3 


1 


.006E- 


-04 


22.4 


1897.0 


1 


.OOOE- 


■04 


40.0 


1754.0 


1 


.OOOE- 


■04 


40.0 


1608.5 


3 


.420E- 


■04 


37.8 


1608.5 


4 


.600E- 


■04 


14.4 


1593.9 


2 


.140E- 


■04 


34.9 


1490.0 


5 


.OOOE- 


■05 


30.0 


1400.0 


3 


.420E- 


■04 


100.0 


1313.0 


1 


.200E- 


■03 


28.0 


1270.0 


1 


.280E- 


■03 


35.0 


1200.0 


3 


.200E- 


■04 


30.0 


1163.1 


8, 


.900E- 


■04 


27.0 


998.0 


1 


.630E- 


■04 


129.1 


940.0 


1 


.600E- 


■04 


13.0 


890.0 


1 


.700E- 


■03 


4.0 


883.0 


1 


.800E- 


■03 


14.1 


836.0 


4 


.280E- 


■04 


14.1 


813.0 


2 


.OOOE- 


■04 


15.0 


800.0 


5 


.300E- 


■04 


70.7 


788.0 


1 


.200E- 


■03 


13.0 


737.0 
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.600E- 


■04 


18.2 


702.0 
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.570E- 


■04 


12.9 


670.0 


8 


.550E- 


■05 


18.3 


635.0 


1 


.710E- 


■04 


18.3 


607.0 


7 


.800E- 


■04 


7.5 


588.0 


8, 


.300E- 


■04 


5.8 


576.0 


4, 


.280E- 


■04 


4.0 


571.0 
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.140E- 


■04 


20.0 


562.0 


8, 


.550E- 


■05 


3.8 


530.0 


1 


.450E- 


■04 


7.0 



1 cm 1^ in normalized units 
(ergcm-2s-lHz-lSr-l) 
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Fig. 1. — The evolution of the spectral energy distribution with age (0 - 10 Myr, every 1 Myr) for a, Z = IZq, logC = 5.0 (log (A^cl) = 5), 
logP/fc = 5.0 Kcm~'^ Hll region-only model (top) and Hll region plus PDR model (bottom). In both diagrams, the final, integrated 
starburst SED is shown as the thick line at high luminosity. 
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X (pm) 

Fig. 2. — PAH emission template. Crosses indicate the observed points, with the soUd Hne indicating our template fit, and the underlying 
curves showing the individual Lorentzian components. 




2 4 6 8 10 

Time (Myr) 

Fig. 3. — The time variation of the incident heating flux (L»/i?gjj) entering the Hll region surrounding a solar metallicity starburst with 
log{Po/fc) = 5. Each curve represents a different C (M^i), and reveals how the C parameterizes different L^/R^^^ at any given time. 
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Fig. 4. — Five SEDs with Solar metallicity, logC = 5.0 and varying pressure (logP/fc = 4, 5, 6, 7 and 8). The model SEDs are almost 
indistinguishable apart from their nebula emission features, such as the [Nell] 12.8 fim emission line. 
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Fig. 5. — Six model SEDs with Solar metallicity, logPo/fc = 5.0 and varying compactness. The compactness parameter decreases from 
logC = 6.5 to logC = 4.0 as the far-IR dust emission feature moves to longer wavelengths. 
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0.1 1.0 tO.O 100.0 

X (/im) 

Fig. 6. — SEDs for a Solar mctallicity starburst with logP/fc = 5, and logC = 5. The top panel has /pdr ranging from 0.0 to 1.0 in 
steps of 0.1. The bottom diagram shows the SEDs using the Tricar parameter introduced in SEDl, with T^icai = 0, 1, 2, 4, 8, 16, and 32 
Myr. It is clear that each of these two formulations of the PDR covering factor problem map to the other. 
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X (pim) 

Fig. 7. — Spectra showing the effect of varying the PDR column depth on a solar metallicity starburst with logC = 5 and logP/fc = 5. 
Shown are a PDR free (Hll region) SED, and PDRs with Hi columns of logAf(HI) = 21.5, 22.0 (standard model), 22.5, and 23.0, as 
labelled. 




Fig. 8. — The UCHII region templates for 0.5, 1.0 and 2.0Zq. These correspond to a Kroupa IMF and a star formation of I.OMq yr~^ 
continued over 1.0 Myr. The small changes in the apparent normalization of these spectra is due to the intrinsic change of stellar luminosity 
with abundance. 
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Fig. 9. — Three SEDs of a solar metallicity starburst with logC = 5, logP/fc = 5 and /pdr = 1-0. Shown is the effect of including the 
IZq UCHII template with /uCHll = 0.0, 0.5, and 1.0, with the mid-IR continuum increasing with increasing /uCHll- 
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Fig. 11. — SED of a starburst with IZq, logC = 5, logP/fc = 5 and /pdr = 10 (lower SED) and SED of the same starburst with our 
IZq "old" stellar SED added (upper SED). We assume a constant SFR of IMQyr^-'^ for the computation of both curves. 




0.1 1.0 10.0 100.0 



Fig. 12. — Attenuation curve from lFischera fc Dopital 1)20051 ) with an Ry = Ay / Eb—v = 3.1. Note that the 2175 A carbon ring tt— orbital 
resonance diffuse absorption feature is weak relative to the extinction curves measured in the Milky Way. Note also the silicate features at 
9.7 and 18 fim. 
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©.1 iM 10.0 100.0 

Fig. 13. — The SED of a model starburst with solar metallicity, lo gC = 5, logP/fc = 5, and /pdr = 0.5 showing the effects of increasing 
attenuation by diffuse dust. Applying the attenuation curve of figure ll2l we show Ay of 0.0, 0.5, 1.0, 2.0, 4.0, 6.0, 10.0, 15.0, 30.0 and 40.0 
as labelled and seen by the decreasing UV and optical flux. Note that the diffuse dust emission associated with the attenuation has not 
been included. 
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Fig. 14. — Cool dust em ission for a Habing radiation field density, Go, of 0.1, 1.0, 10.0, 100.0 times the Local Interstellar Radiation Field 
Density of IHabin^ I I1968I) . The strength of the radiation field increases from bottom to top. 
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Fig. 15. — Five SEDs demonstrating the effect of adding diffuse dust emission at 10% of the Starburst luminosity to our fiducial model 
Starburst (IZq, logC = 5, logP/fc = 5, and /pdr = 0.5). The SEDs show the original SED (thick line) and those with an added diffuse 
dust continuum heated by a young stellar continuum 100.0, 10.0, 1.0 and 0.1 times that of the local ISRF. 
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Fig. 17. — The fit to the starburst galaxy NGC 6240. The errors in the observed fluxes are similar in magnitude to the height of the 
blue squares. The parameters of the fit are given on the label. Note that, although here Ay is estimated at 1.9 mag, this applies to the 
older stars. For the younger stars in the Hll regions, Ay ~ 4.4 mag in this model. 
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Fig. 18. — The sensitivity of fit to tiie starburst galaxy NGC 6240 to the various parameters. As can be seen, for this galaxy, logC is 
constrained within ±0.25 dex. Ay to ±0.3 mag, and f^id to within 20%. The fraction of UCHII regions, /uCHll is not very well constrained 
in this galaxy, owing to its high compactness, but would be much better constrained in galaxies with logC < 5. 



32 




33 



44,0 



43.0 



42,0 



41,0 



40-0 



39.0 



n-r-r 



T 



T 



Arp2£0 

I I I I I r I I I I I { 



- Arp220 : SFR=31 5; T^OA; log C=5-0; A(V)(SB) = 2Q 
A(V)=2.2; F{old)=1 .0 ; F(UCHI1)= 0.7 




log ?^ (jim) 



Fig. 20. — The double foreground dust screen fit to the starburst galaxy Arp 220. The screen around the starburst nucleus has 
Av = 20 mag and the diffuse older star component has an attenuation of Ay = 2.2. The fit around the silicate absorption features is 
improved in this more complex model. 
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